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Using a Coarse-grained Molecular Dynamics (CMD) approach we study the apparent nonlinear 
dynamics of water molecules filling/emptying carbon nanotubes as a function of system parameters. 
Different levels of the pore hydrophobicity give rise to tubes that are empty, water- filled, or fluctu- 
ate between these two long-lived metastable states. The corresponding coarse-grained free energy 
surfaces and their hysteretic parameter dependence are explored by linking MD to continuum fixed 
point and bifurcation algorithms. The results are validated through equilibrium MD simulations. 
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Molecular dynamics (MD) simulations on classical or 
quantum energy surfaces provide an atomically detailed 
description of complex molecular processes such as pro- 
tein folding or materials fracture. However, such simu- 
lations are inherently limited by the requirement to in- 
tegrate accurately even the fastest molecular processes 
of bond vibrations and atomic collisions. Remarkable 
progress has been made recently in addressing the time 
scale limitations of MD (see, e.g., |2l, La Li)- To inte- 
grate out the fast molecular motions and explore the slow 
dynamics, we are developing a "coarse molecular dynam- 
ics" (CMD) approach 0, Q. In principle, the projection 
operator formalism Q connects the microscopic dynam- 
ics to such slow motions. However, exact analytical con- 
struction of the corresponding noise and memory terms 
is usually intractable. In CMD, we circumvent this chal- 
lenge by estimating on the fly both the thermodynamic 
driving forces for slow motions and their dynamic prop- 
erties. Information about the slow coarse dynamics is 
extracted from the projected motions of many, appropri- 
ately initialized, but otherwise independent and unbiased 
replica simulations. 

In this letter, we extend the CMD approach to ex- 
plore directly the physical parameter space of a com- 
plex molecular system. This is accomplished by link- 
ing short bursts of MD simulations with continuum al- 
gorithms to simultaneously search phase and parameter 
space. The general formalism is illustrated by study- 
ing the parametrically controlled equilibrium between the 
water-filled and empty state of a short and narrow car- 
bon nanotube (CNT) dissolved in water (Figure [lj. The 
CNT in water serves as a paradigm of a molecular sys- 
tem with nontrivial dynamic and thermodynamic behav- 
ior under parametric control. Conventional MD simula- 
tions showed that for certain interaction strengths be- 
tween the carbon atoms and water, the molecular system 
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FIG. 1: Schematic of the CNT- water system showing the 
"vestibules" at the openings of the CNT and corresponding 
dimensions. The CNT is represented as a cylinder surrounded 
by water molecules. 



exhibits "bi-phasic" character, with the tube fluctuating 
on a nanosecond time-scale between a water-filled and an 
empty state. Such controlled fluctuations in the water oc- 
cupancy have been suggested as a mechanism to regulate 
water, proton,and ion transport through biological pro- 
tein pores @, • The computational prediction of filled 
tubes at strong water-tube interactions has since been 
confirmed experimentally for somewhat wider tubes [HI . 
Important here is that, despite the molecular complexity, 
the equilibrium thermodynamic calculations required to 
validate CMD are feasible. 

We perform classical MD simulations of a (6,6)-type 
nanotube of ~13.5 A length and 8.1 A diameter in a 
box of N wat = 1034 TIP3P water molecules [ll| under 
periodic boundary conditions with Ewald electrostatics 
and a time step of 0.002 ps (additional simulation de- 
tails as in Refs. || 0). We use CMD to explore the 
thermodynamics and kinetics of tube filling as a func- 
tion of the "hydrophobicity parameter" A scaling the 
attractive r~ 6 carbon-water Lennard- Jones interactions: 
mc-oW = Ar~ 12 - XCr~ 6 where A = 696,790.7 kcal 
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mol -1 A 12 and C = 564.5036 kcal mol -1 A 6 , respectively. 
In physical experiments, the effect of A on the water- 
affinity for the CNT pore can be mimicked by changing 
quantities like water pressure or osmolality. To mon- 
itor filling and emptying, we use the water occupancy 
as a coarse observation variable, N = 5^^ at Q(ri, Zi), 
where Ti and zi are the water radial and axial positions in 
the cylindrical coordinate system defined by the instanta- 
neous position and orientation of the freely moving CNT 
(other simulation details as in 8]). The weight function 
is given by 9(r, z) = exp[- (r / R cy if o - (2z/L cy i) 6 } where 
L cy i = 17.5 A and R cy i = 4.05 A are the length and 
radius, respectively, of a cylindrical region that extends 
somewhat beyond the ~13.5 A long CNT to include wa- 
ter molecules near the openings (see FigureGJ. By includ- 
ing these "vestibules" of the CNT, we resolve the gradual 
entry, or exit, of a water molecule without compromising 
the information about the state of nanotube filling. We 
found earlier that the tube appears bi-stable for A ~ 0.8 
Et : it is predominantly empty below A « 0.7, and 
filled above A « 0.9. 

Bi-stability can be succinctly summarized in a bifur- 
cation diagram that reports local maxima of the equilib- 
rium density as a function of A and - in the bi-stable 
regime - also the intervening saddles. In what follows, 
we will use short replica MD simulations to construct 
a kinetically-based coarse bifurcation diagram and com- 
pare it to thermodynamics. We will report the fixed 
points of a timestepper constructed by CMD through 
the following steps. In the "lifting" step, we create rep- 
resentative initial phase points of the full molecular sys- 
tem. For a given value of the hydrophobicity param- 
eter A, we harmonically bias the occupancy N during 
a short, 15-ps MD run, toward a target No by adding 
Unas = k(N - N ) 2 / 2 to the Hamiltonian (k = 50 kcal 
mol -1 for N > 1; k = 80 kcal mol" 1 for N < 1). 
Five different conformations near each target No value 
are saved during the last 5 ps of the constrained run, 
and 10 sets each of random Maxwell-Boltzmann veloci- 
ties are assigned to create 50 starting configurations. In 
the "evolve" step, an unconstrained MD run of r = 1 ps 
duration is performed for each structure. Finally, in the 
"restrict" step, the resulting trajectories are projected 
onto the coarse variable N. In particular, the average 
occupancy N at t he end of th e short runs defines our 
CMD timestepper: N(r; N ; A) = $ r (iVo; A). 

The fixed points for given values of A satisfy 



f(N; X) = N — <& r (N; A) = 0. 



(1) 



They correspond to metastable empty or full, as well as to 
unstable, partially filled "transition states" . For small r 
the above equation approximates the differential (steady 
state) form d$t(No; X)/dt = 0; all our coarse-grained 
computations can be equivalently performed through this 
differential form. Fixed points for given A are computed 
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FIG. 2: Coarse bifurcation diagram for the CNT- water sys- 
tem. Solid filled squares (■) correspond to fully filled (N > 5) 

and empty states (N < 1), respectively. Open circles (O) 
correspond to partially filled states, while the filled diamonds 
(♦) correspond to the turning points of the coarse bifurca- 
tion diagram. Vertical dashed lines indicate the upper (A) 
and lower (F) turning points, and A = 0.85 (points C-D- 
E), respectively. Representative structures for points A-F are 
indicated by arrows (carbon: gray; oxygen: red; hydrogen: 
white). The coarse effective free energy surfaces for A = 0.7, 
0.85, and 0.95, respectively, are shown below the main figure. 
The inset in the top right corner of the main figure shows 
the effective bifurcation diagram obtained through histogram 
reweighting. 

by solving Eq. (Q through a Newton- Raphson type iter- 
ation with an appropriate initial guess for N: 



Ni=N 



f(Np;X) 
df/dNo 



(2) 



In Eq. (0), df/dNo is estimated by computing N — 
& T (N;X) for nearby initial points N. Fixed points con- 
verged to within 0(1O -3 ). The results (fixed points as a 
function of A) are plotted in Figure [2j 

At the upper (lower) turning points we lose the filled 
(empty) fixed points. In principle, turning points can be 
found by computing the entire diagram (using pseudo- 
arclength continuation Q); however, numerical bifurca- 
tion theory provides more efficient algorithms that solve 
directly for these points (and thus the boundaries of ap- 
parent hysteresis in parameter space). In compact no- 
tation, X = [A/ - , A] T is a vector describing the position 
in the bifurcation diagram; the turning points are roots 
of the vector- valued function F(X) = [N - & t (N; A), 1 - 
d$ t (N; X)/dN] T = 0. To solve these two coupled equa- 
tions for the critical N and A simultaneously, we use a 
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2D recursive Newton-Raphson iteration procedure: 
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dF 
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The derivatives of the Newton-type iteration, 
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are estimated numerically from replica simulations ini- 
tiated on a centered "stencil" (i.e., grid) with \AN\ = 
0.1 and |AA| = 0.025 around X = [N ,X } T . The 
"lift-evolve-restrict" procedure is implemented and the 
timesteppers, <3> T (iV;A), are computed for each node of 
the grid. Starting from a good initial guess, the proce- 
dure converges after a few (~5) iteration steps: the up- 
dates for [TV, A] T and residuals for F are of order O(10~ 3 ) 
and 0(1O -2 ), respectively. 

Figure[2]shows the resulting coarse bifurcation diagram 
for the CNT- water system. The two metastable branches 
of the S-shaped diagram correspond to empty (N « 1) 
and filled states (N w 7). The lower and upper turning 
points are at A ~ 0.7 and 0.95, respectively. 

So far, we have taken an entirely dynamic perspective 
in our computations. To estimate the underlying thermo- 
dynamics, knowing that N is a good reaction coordinate 
[l3|, we assume diffusive dynamics along N in the form 
of an effective Fokker-Planck (FP) equation: 



dP(N,t;\) 
dt 



-^(JV;A) + ^23(JV;A)]P(JV,t;A). 

(5) 

The coarse timestepper calculations are precisely what is 
needed to parametrize this FP: 



v(N;X) . ^V^; 2D(N; A) = 



dt 



dt 



(6) 

For short simulation times r, the coarse timestepper (the 
average of the results of the replica simulations) is used 
to estimate the drift v(N; A); from the time-dependence 
of the variance a 2 [N(t; No; A)] of TV, we estimate the dif- 
fusion coefficient D(N] A). In this notation, the effective 
free energy G(AT; A) for the problem is 



0G(N; A) 



Jo 



v(N';\) 
D(N';\) 



dN' + In D(N; A) + const. 

(7) 

Maxima of the equilibrium density exp(— /3G(N; A)) can 
be found as zeroes of v(N;\) — D'(N;\) [or, ignor- 
ing a weak N dependence of D, by the zeroes of 
v(N; A)]. Clearly, the fixed points of our coarse timestep- 
per are approximations of the free energy extrema, where 
metastable and unstable fixed points correspond to well 
bottoms and intervening saddles, respectively. 

In the inset of Figure [2j we show for comparison the 
thermodynamic bifurcation diagram, indicating the min- 
ima and saddles of the equilibrium free energy surfaces 



G(N; A). This diagram is estimated from three long (132 
ns total) equilibrium MD simulations, one in the pre- 
dominantly filled regime (A = 1), one in the predom- 
inantly empty regime (A ~ 0.75), and one in between 
(A = 0.785). By combining the 3-dimensional histograms 
of occupancy fluctuations and water- CNT interaction en- 
ergies using a weighted histogram method 0, we 
obtain the free energy surface G(N;X) at different A val- 
ues. To avoid difficulties with insufficient sampling of the 
observation variable in the equilibrium runs, we coarse- 
grained TV to the nearest integer numbers. 

Overall, we find excellent agreement between the 
coarse bifurcation diagram from multiple short (1-ps) 
CMD runs, and that from long equilibrium runs. In 
particular, the location of the metastable and unstable 
branches is nicely reproduced by the equilibrium results. 
The main difference is that the lower turning point is 
at a somewhat lower interaction strength A « 0.66 com- 
pared to 0.7. Such small differences can be due to sev- 
eral factors: (1) The equilibrium-predicted lower turn- 
ing point lies outside the A-range probed by MD and 
had to be extrapolated. (2) The equilibrium computa- 
tions were coarse-grained to integer N values. (3) The 
N- dependence of the effective diffusion coefficient was ne- 
glected in constructing the kinetic bifurcation diagram; 
and finally (4) CMD simulations at stationarity probe 
the dynamic stability of an apparent fixed point over a 
given time horizon r. The location of fixed as well as 
turning points from CMD will depend somewhat on r. 
Here, the short runs are about 1 ps long, one order of 
magnitude shorter than the duration of a typical wa- 
ter uptake or release event N ^ N + 1 assuming 
Markov-chain dynamics. For parameter values close to 
the turning points, as the free energy wells become more 
shallow, typical transition times N ^ 7V+1 become com- 
parable to any fixed r. As a result, one would expect the 
CMD turning points to move "inwards," and the appar- 
ent hysteresis range to shrink [lf| . 

Our CMD calculations have so far been local in N-X 
space. As an example of a parametric search involving 
global properties, we next locate the critical hydropho- 
bicity parameter value at which the filled and empty 
states have (approximately) equal populations. We for- 
mulate this as a fixed point problem that we solve it- 
eratively. Given an initial guess of A we first solve for 
the two metastable fixed points [Af emp t y (A) and iVf u n(A)]. 
We then consider any reasonable path between the two 
(parametrized by the coarse observable A 7 "), discretize it, 
and estimate (through locally initialized and executed 
MD experiments and quadrature) the state dependent 
drift and diffusion coefficient in the FP equation, Eq. (0). 
Finally, using Eq. 0, we compute G(N; A). To solve the 
equation AG(A) = G(A^ ull ; A) - G(Af empty ; A) = for the 
critical A we use a contraction mapping with numerically 
estimated derivatives. G(N; A) was computed as follows: 
From different initial configurations for < A/^ < 8, we 
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FIG. 3: Coarse free energy surface as a function of N 
for A = 0.785 where the empty and filled wells have equal 
depth. Results from equilibrium MD (dashed) and CMD 
[from Eq. Q; solid] are shown. The effective free energy pre- 
dicted by CMD agrees well with that of 19-ns long equilibrium 
MD simulations. 



run 50 short (r = 1 ps) replica simulations. By robust lin- 
ear fits to the evolution of N(t; Ni\ A) and cr 2 [N(t; iV*; A)], 
we estimate v(N) and D(N) and compute (3G(N) from 
Eq. by integration through a smooth spline approxi- 
mation to the data. As illustrated in Figure the criti- 
cal A value is approximately 0.785, and this is confirmed 
by equilibrium MD simulations. In particular, the MD 
and CMD free energy surfaces at this critical parame- 
ter value are in excellent agreement. We note that for 
equal populations in the two wells, the free-energy gradi- 
ents estimated by CMD are relatively small, resulting in 
particularly reliable free energy surfaces, with larger de- 
viations expected in regimes where one or the other well 
dominates. Finally, we note that CMD in combination 
with Kramers theory for diffusive barrier crossing, pro- 
duces accurate estimates of the rate coefficients for escape 
from the metastable filled and empty states, when com- 
pared with equilibrium MD [~1/(140 ps) vs. ~1/(190 
ps) from MD at the transition midpoint, A = 0.785]. 

We illustrated the extension of the CMD approach 
to the exploration of parameter space of molecular sys- 
tems. Our computational methods are general and can be 
used to study metastability boundaries in other systems 
(e. g., k inetic Monte Carlo 0, ^3 or Brownian Dynam- 
ics Linking continuum nonlinear dynamics compu- 
tation (fixed point, continuation, integration or bifurca- 
tion algorithms) with MD simulations has the potential 
to accelerate parametric exploration of dynamic as well 



as thermodynamic features. The knowledge of good re- 
action coordinates is crucial for the success of this ap- 
proach; determining such coordinates based on data is 
the subject of intense current research [3. Hfj Eol . 
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